#include <iostream>
#include <vector>
#include <numeric>
#include <omp.h>
#include <random>

int main(int argc, char** argv) {
  float p=atof(argv[1]);
  int N=atoi(argv[2]);
  std::vector<int> count;
  std::cout << "probability " << p << std::endl;

  int T = omp_get_max_threads();
  std::cout << "Threads " << T << std::endl;

  {
	std::random_device rd;
	std::mt19937 gen(rd());
	std::uniform_real_distribution<> dis(0, 1);

	int i=0;
	while (i<N) {
	  int k=0;
	  int v=(dis(gen)>=p)?0:1;
	  i++;
	  while (v==1) {
		k++;
		v=(dis(gen)>=p)?0:1;
		i++;
	  }
	  if (k) {
		{
		  count.push_back(k);
		}
	  }
	}
  }
  int sum = std::accumulate(count.begin(), count.end(),0);
  std::cout  << "count = " << count.size() 
			 << " Average " 
			 << static_cast<double>(sum)/static_cast<double>(count.size()) 
			 <<  " Est. " << 1./(1.-p) << std::endl;

  return 0;
}
				 
